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High Resolution Methods 
for Time Dependent Problems 
with Piecewise Smooth Solutions 



Eitan Tadmor* 

Abstract 

A trademark of nonlinear, time-dependent, convection-dominated prob- 
lems is the spontaneous formation of non-smooth macro-scale features, like 
shock discontinuities and non-differentiable kinks, which pose a challenge for 
high-resolution computations. We overview recent developments of modern 
computational methods for the approximate solution of such problems. In 
these computations, one seeks piecewise smooth solutions which are realized 
by finite dimensional projections. Computational methods in this context can 
be classified into two main categories, of local and global methods. Local 
methods are expressed in terms of point-values ( — Hamilton-Jacobi equa- 
tions), cell averages ( — nonlinear conservation laws), or higher localized mo- 
ments. Global methods are expressed in terms of global basis functions. 

High resolution central schemes will be discussed as a prototype example 
for local methods. The family of central schemes offers high-resolution "black- 
box-solvers" to an impressive range of such nonlinear problems. The main in- 
gredients here are detection of spurious extreme values, non-oscillatory recon- 
struction in the directions of smoothness, numerical dissipation and quadra- 
ture rules. Adaptive spectral viscosity will be discussed as an example for 
high-resolution global methods. The main ingredients here are detection of 
edges in spectral data, separation of scales, adaptive reconstruction, and spec- 
tral viscosity. 
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1. Introduction — piecewise smoothness 
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A trademark of nonlinear time-dependent convection-dominated problems is 
the spontaneous formation of non-smooth macro-scale features which challenge 
high-resolution computations. A prototype example is the formation of shock dis- 
continuities in nonlinear conservation laws, 

d 

—u(x,t)+V x -f(u(x,t))=0, u:= (u u ...,u m ) T . (1.1) 

It is well known, e.g., ||], that solutions of Hl.lfl cease to be continuous, and 
should be interpreted in a weak sense with the derivatives on the left as Radon 
measures. This requires clarification. If u(x,t) and v(x,t) are two admissible solu- 
tions of l|l.lfl then the following stability estimate is sought (here and below a, f3, ... 
stand for different generic constants), 

\H;t)-v(;t)\\<a t \\u(;0)-v(;0)\\. (1.2) 

Such estimates with different norms, || • ||, are playing a key role in the linear 
setting — both in theory and computations. For linear hyperbolic systems, for 
example, (|1.2|) is responsible for the usual L 2 -stability theory, while the stability of 
parabolic systems is often measured by the L°°-norm, consult . But for nonlinear 
conservation laws, (|1.2|) fails for any L p -norm with p > 1. Indeed, comparing u(x, t) 
with any fixed translation of it, v[x, t) :— u(x + h,t), the L p version of i|1.2(l implies 

\\A +h u(-,t)\\ LP{Rd) < a t \\A +h u(-,0)\\ LP{Rd) , A +h u(-,t) := u(- +h,t)- u(-,t). 

For smooth initial data, however, the bound on the right yields ||A_|_fe'u(-, i)||z,p < 
Q!f|/i|, which in turn, for p > 1, would lead to the contradiction that u(-,t) must 
remain continuous. Therefore, conservation laws cannot satisfy the L p -stability es- 
timate <|1.2I) after their finite breakdown time, except for the case p = 1. The latter 
leads to Bounded Variation (BV) solutions, ||u(-,t)||sv := sup ||A+/ l 'u(',t)||ii/|/i| < 
at < oo, whose derivatives are interpreted as the Radon measures mentioned above. 
BV serves as the standard regularity space for admissible solutions of Ijl.lfl . A com- 
plete BV theory for scalar conservation laws, m = 1, was developed Kruzkov. Fun- 
damental results on BV solutions of one-dimensional systems, d = 1, were obtained 
by P. Lax, J. Glimm, and others. Consult for recent developments. Relatively 
little is known for general (m — 1) X (d — 1) > 0, but cf., |12| . 

We argue that the space of BV functions is still too large to describe the ap- 
proximate solutions of (|l.l(l encountered in computations. Indeed, in such compu- 
tations one does not 'faithfully' realize arbitrary BV functions but rather, piecewise 
smooth solutions. We demonstrate this point in the context of scalar approximate 
solutions, v h (x,t), depending on a small computational scale h ~ 1/N. A typical 
error estimate for such approximations reads, 0j 

\\v h (;t) - u(;t)\\Ll e (R) < 0) - u(-, 0) H^JR) + « t ^ /2 . (1.3) 

The convergence rate of order 1/2 is a well understood linear phenomena, which is 
observed in computations 1 . The situation in the nonlinear case is different. The 

Bernstein polynomials, B^(u), provide a classical example of first-order monotone approxi- 
mation with iVerror of order (||u||fly /.AQ 1 / 2 - The general linear setting is similar, with improved 
rate ~ /i r /( r + 1 ) for r-order schemes. 
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optimal convergence rate for arbitrary BV initial data is still of order one-half, 
|15|. but actual computations exhibit higher-order convergence rate. The apparent 
difference between theory and computations is resolved once we take into account 
piecewise smoothness. We can quantify piecewise smoothness in the simple scalar 
convex case, where the number of shock discontinuities of u(-,t) is bounded by the 
finitely many inflection points of the initial data, u(x, 0). In this case, the singular 
support of it(-, t) consists of finitely many points where S{t) = {x \ d x u(x, t) J, — oo}. 
Moreover, the solution in between those point discontinuities is as smooth as the 
initial data permit, |23i namely 

sup \&Pu(x,t)\ < e pLT sup \d%u(x,0)\+Const L , S L {t) := {x \ d x u(x,t) > -L}. 

x<ES L (t) x£S L {0) 

If we let d(x,t) := dist(x,S(t)) denote the distance to S(t), then according to |19| . 
the following pointwise error estimate holds, \v (x,t) — u(x,t)\ < ath/d(x,t), and 
integration yields the first-order convergence rate 

H^C-.t) - t*(-,*)|| i i )cC R ) < at/i|log(/i)|. (1.4) 

There is no contradiction between the optimality of 1)1. 3[) and 11. 4|) . The former ap- 
plies to arbitrary BV data, while the latter is restricted to piecewise smooth data and 
it is the one encountered in actual computations. The general situation is of course, 
more complicated, with a host of macro-scale features which separate between re- 
gions of smoothness. Retaining the invariant properties of piecewise smoothness in 
general problems is a considerable challenge for high-resolution methods. 

2. A sense of direction 

A computed approximation is a finite dimensional realization of an underlying 
solution which, as we argue above, is viewed as a piecewise smooth solution. To 
achieve higher accuracy, one should extract more information from the smooth parts 
of the solution. Macro-scale features of non-smoothness like shock discontinuities, 
are identified here as barriers for propagation of smoothness, and stencils which 
discretize 1)1.111 while crossing discontinuities are excluded because of spurious Gibbs' 
oscillations. A high resolution scheme should sense the direction of smoothness. 

Another sense of directions is dictated by the propagation of information gov- 
erned by convective equations. Discretizations of such equations fall into one of 
two, possibly overlapping categories. One category of so-called upwind schemes 
consists of stencils which are fully aligned with the local direction of propagating 
waves. Another category of so-called central schemes consists of two-sided stencils, 
tracing both right-going and left-going waves. A third possibility of stencils which 
discretize 'against the wind' is excluded because of their inherent instability, 

|14| . A stable scheme should sense the direction of propagation. 

At this stage, high resolution stable schemes should compromise between two 
different sets of directions, where propagation and smoothness might disagree. This 
require essentially nonlinear schemes, with stencils which adapt their sense of direc- 
tion according to the computed data. We shall elaborate the details in the context 
of high-resolution central scheme. 
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3. Central schemes 

We start with a quotation from ^] §12.15], stating "In 1959, Godunov de- 
scribed an ingenious method for one-dimensional problems with shocks" . Godunov 
scheme is in the crossroads between the three major types of local discretizations, 
namely, finite-difference, finite-volume and finite-element methods. The ingenuity 
of Godunov's approach, in our view, lies with the evolution of a globally defined 
approximate solution, v (x,t n ), replacing the prevailing approach at that time of 
an approximate solution which is realized by its discrete gridvalues, v v {t n ). This 
enables us to pre-process, to evolve and to post-process a globally defined approx- 
imation, v h (x,t n ). The main issue is how to 'manipulate' such piecewise smooth 
approximations while preserving the desired non-oscillatory invariants. 

Godunov scheme was originally formulated in the context of nonlinear con- 
servation laws, where an approximate solution is realized in terms of a first-order 
accurate, piecewise-constant approximation 

v h (x,t n ) := A h v(x,t n ) :=y^v v {t n )l Iv {x), v v {t n ) := ±- [ v(y,t n )dy. 

The cell averages, v v , are evaluated over the equi-spaced cells, I v := {x \ \x — x v \ < 
h/2} of uniform width h = Ax. More accurate Godunov- type schemes were devised 
using higher-order piecewise-polynomial projections. In the case of one-dimensional 
equi-spaced grid, such projections take the form 

v h v(x) = ^Tp v (x)ii v (x), p v (x) = v v + v v{—-f^~) + \ v "(c~jr^~) +•••• 

V 

Here, one pre-process the first-order cell averages in order to reconstruct accurate 
pointvalues, v v , and say, couple of numerical derivatives vj /h,vj' /h 2 , while the 
original cell averages, {v u }, should be preserved, AhPhV h = Ai x v h . The main 
issue is extracting information in the direction of smoothness. For a prototype 
example, let A+v v and A-V v denote the usual forward and backward differences, 
A±v v := ±{v v ±i — v v ). Starting with the given cell averages, {v„}, we set v v = v v , 
and compute 

/ /a - a - \ / \ sgn(zi) + sgn(z 2 ) . ,, , , n /olN 
v u = mm(A + Vv, A_v„), mm(z 1 , z 2 ) ■= mui{|^i|, |z 2 |}. (3.1) 

The resulting piecewise-linear approximation is a second-order accurate, Total Vari- 
ation Diminishing (TVD) projection, \\T'hV h {x)\\BV < ||w' l (a;)||By. This recipe of 
so-called minmod numerical derivative, l|3.1|l . is a representative for a large library 
of non-oscillatory, high-resolution limiters. Such limiters dictate discrete stencils in 
the direction of smoothness and hence, are inherently nonlinear. Similarly, nonlinear 
adaptive stencils are used in conjunction with higher-order methods. A description 
of the pioneering contributions in this direction by Boris & Book, A. Harten, B. 
van-Leer and P. Roe can be found in [TO] . The advantage of dealing with globally 
defined approximations is the ability to pre-process, to post-process and in partic- 
ular, to evolve such approximations. Let u{x, t) = u h (x,t) be the exact solution 
of (|1.1|) subject to u h (x,t n ) — VhV h (x,t n ). The exact solution lies of course out- 
side the finite computational space, but it could be realized in terms of its exact 
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cell averages, v h (x, t n+1 ) = ^ v v v (t n+1 )li i/ (x). Averaging is viewed here a simple 
post-processing. Two prototype examples are in order. 

Integration of l|l.l|l over control volume l v x [t n ,t n+1 ], forms a local stencil which 
balances between the new averages, {v v (t n+1 )}, the old ones, {v v+ k(t n )} , and the 
fluxes across the interfaces along x„±i/2 x [t n ,t n+1 ]. In this case, the solution along 
these discontinuous interfaces is resolved in terms of Riemann solvers. Since one 
employs here an exact evolution, the resulting Godunov-type schemes are upwind 
scheme. The original Godunov scheme based on piecewise constant projection is 
the forerunner of all upwind schemes. As an alternative approach, one can realize 
the solution u(x,t n+1 ), in terms of its exact staggered averages, {v u+ i/ 2 (t n+1 )}. 
Integration of (|l.l|l over the control volume /„+i/2 x [t n ,t n+1 ] subject to piecewise 
quadratic data given at t = t n , u(x,t n ) — VhV h (x,t n ), yields 

+ ~ (vj(n v' v+ i(n) + H [Kil 12 k +1/2 ) , (3.2) 

where jt , ™+ 1 / 2 stands for the averag ed flux, Fl l+1/2 = r t T f[u{x u ,T)dr/M, 
Thanks to the staggering of the grids, one encounters smooth interfaces x v x 
[i n ,i n+1 ], and the intricate (approximate) Riemann solvers are replaced by sim- 
pler quadrature rules. For second-order accuracy, for example, we augment (|3.2fl 
with the mid-point quadrature 

^ n+1/2 = /M*" +1/2 )), v v (t n+1/2 ) = Mn-^fMt n ))'- (3-3) 

z/\x 

Here, the prime on the right is understood in the usual sense of numerical dif- 
ferentiation of a gridfunction - in this case the flux {f(v u (t n ))}i,- The resulting 
second-order central scheme (|3.3H . (|3.2[1 was introduced in ^3]. It amounts to a sim- 
ple predictor-corrector, non-oscillatory high-resolution Godunov-type scheme. For 
systems, one implements numerical differentiation for each component separately. 
Discontinuous edges are detected wherever cell-averages form new extreme values, 
so that vj(t n ) and f(v v (t n ))' vanish, and (|3.3|l . (|3.2() is reduced to the forerunner of 
all central schemes — the celebrated first-order Lax-Friedrichs scheme, This 
first-order stencil is localized to the neighborhood of discontinuities, and by assump- 
tion, there are finitely many them. In between those discontinuities, differentiation 
in the direction of smoothness restores second-order accuracy. This retains the 
overall high-resolution of the scheme. Consult Figure [I] for example. 

Similarly, higher-order quadrature rules can be used in connection with higher- 
order projections, H2, |llj . A third-order simulation is presented in Figure|21 Finite- 
volume and finite-element extensions in several space dimensions are realized over 
general, possibly unstructured control volumes, f2„ x [t n , t n+1 ], which are adapted to 
handle general geometries. Central schemes for 2D Cartesian grids were introduced 
in and extended to unstructured grids in pQ. A similar framework based on 
triangulated grids for high-resolution central approximations of Hamilton- Jacobi 
equations was described in [3] and the references therein. 

Central schemes enjoy the advantage of simplicity - they are free of (approx- 
imate) Riemann solvers, they do not require dimensional splitting, and they apply 
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to arbitrary flux functions 2 without specific references to eigen-decompositions, Ja- 
cobians etc. In this context, central schemes offer a "black-box solvers" for an 
impressive variety of convection-dominated problems. At the same time, the cen- 
tral framework maintain high-resolution by pre -and post-processing in the direc- 
tion of smoothness. References to diverse applications such as simulations of semi- 
conductors models, relaxation problems, geometrical optics and multiphase compu- 
tations, incompressible flows, polydisperse suspensions, granular avalanches MHD 
equations and more can be found at jlfij . 




Figure 1 : Second-order central scheme 
simulation of semiconductor device 
governed by ID Euler-Poisson equa- 
tions. Electron velocity in 10 7 cm/s 
with N = 400 cells. 



Figure 2: Third-order central scheme 
simulation of ID MHD Riemann prob- 
lem, with N = 400 cells. The y- 
magnetic field at t = 0.2. 



The numerical viscosity present in central schemes is of order — ■ It is 

suitable for the convective regime where At ~ Ax, but it is excessive when a small 
time step is enforced, e.g., due to the presence of diffusion terms. To overcome this 
difficulty, a new family of central schemes was introduced in 7 and was further 
refined in H]. Here, the previous staggered control volumes, I u +i/2 x [* n j*" +1 ]j 
is replaced by the smaller — and hence less dissipative, J v +i/2 X [t n ,t n+1 ], where 
Ju+i/2 '■= x v+i/2 + At x [a~,a + ] encloses the maximal cone of propagation, = 
a v+\/2 ~ (min) k^k (/«) • The ^ act that tne staggered grids are O(At) away from 
each other, yields central stencils with numerical viscosity of order 0(Ax 2r ^ 1 ). 
Being independent of At enables us to pass to the limit At J, 0. The resulting semi- 
discrete high-resolution central scheme reads v v (t) — —{fv+xnif) ~ fu-i/2(t))/ Ax, 
with a numerical flux, /„+i/2, expressed in terms of the reconstructed pointvalues, 

V ± = V - + l/2 = T D hV h (x l , + 1 /2±,t), 

U+i/2{t):= = +a^a — -. 3.4) 



2 An instructive example is provided by gasdynamics equations with tabulated equations of 
state. 
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Figure 3: Convectio n-diffusi on eq. Figure 4: Third-order central scheme 
u t + (u 2 ) x /2 = (u x /y/l + ul) x sim- for 2D Riemann problem. Density con- 
ulated by , (£0} , with 400 cells. tour lines with 400 x 400 cells. 



4. Adaptive spectral viscosity methods 

Godunov-type methods are based on zeroth-order moments of In each 

time step, one evolves one piece of information per spatial cell — the cell aver- 
age. Higher accuracy is restored by numerical differentiation in the direction of 
smoothness. An alternative approach is to compute higher-order moments, where 
the cell averages, v v , and say, couple of numerical derivatives, vj, v v " are evolved in 
time. Prototype examples include discontinuous Galerkin and streamline-diffusion 
methods, where several local moments per computational cell are evolved in time, 
consult 0] and the references therein. As the number of projected moments in- 
crease, so does the size of the computational stencil. At the limit, one arrives at 
spectral methods based on global projections, 

v N (x,t) = Vpfv(x,t) := ^ v k (t)4> k (x), Vk •= (v(-,t),(p k (-))- 

\k\<N 

Here, 4> k {x) are global basis functions, 4>k — e lkx , Tk(x), ... and v k are the mo- 
ments induced by the appropriately weighted, possibly discrete inner-product (•, •). 
Such global projections enjoy superior resolution — the error \\v(-) — Vnv(')\\ de- 
cays as fast as the global smoothness of v(-) permits. With piecewise smooth so- 
lutions, however, we encounter first-order spurious Gibbs' oscillations throughout 
the computational domain. As before, we should be able to pre- and post-process 
piecewise smooth projections, recovering their high accuracy in the direction of 
smoothness. To this end, local limiters are replaced by global concentration ker- 
nels of the form, K v N v N {x) = J2\k\<N 'Miv) Vk((/>k( x ))x, where r?(-) is an 
arbitrary unit mass C°° [0, 1] function at our disposal. Detection of edges is fa- 
cilitated by separation between the 0(1) scale in the neighborhoods of edges and 
0{h r ) scales in regions of smoothness, Knvn(x) ~ [v(x)) + 0{Nd(x))~ r . Here, 
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[v(x)] denotes the amplitude of the jump discontinuity at x (- with vanishing am- 
plitudes signaling smoothness), and d(x) is the distance to the singular support 
of v(-). Two prototype examples in the 27r-periodic setup are in order. With 
= 1 one recovers a first-order concentration kernel due to Fejer. In 5 we 
introduced the concentration kernel, f? e:cp (£) :~ exp{(3/S,(l — £)}, with exponen- 
tially fast decay into regions of smoothness. Performing the minmod limitation, 
<|3.1|l . mm{KjfVN(x), K^ p vn(x)}, yields an adaptive, essentially non-oscillatory 
edge detector with enhanced separation of scales near jump discontinuities. Once 
macro-scale features of non-smoothness were located, we turn to reconstruct the 
information in the direction of smoothness. This could be carried out either in the 
physical space using adaptive mollifiers, \&e, p , or by nonlinear adaptive filters, a p . 
In the 2n periodic case, for example, we set * Vnv(x) := {^>e, P {x — -),Vnv(-)) 
where \& is expressed in terms of the Dirichlet kernel, D p (y) :— s ™ ' 

*oAv) = \p(S) d A\ ) ' 9 = My) := e ^ 2/fe2_1)l [-M](y)- (4-i) 

Mollifiers encountered in applications maintain their finite accuracy by localization, 
letting 9 { 0. Here, however, we seek superior accuracy by the process of cancellation 

with increasing p ]. To guarantee that the reconstruction is supported in the 

direction of smoothness, we maximize 9 in terms of the distance function, 9(x) = 

d(x), so that we avoid crossing discontinuous barriers. Superior accuracy is achieved 

by the adaptive choice p :~ d(x)N, which yields, [23! , 

\*{d( x ),dWN} * 'Pnv(x) - v(x)\ < Const.(d(x)N) r e^V^. 

The remarkable exponential recovery is due to the Gevrey regularity of pp € Qi and 
Figure [S] demonstrates such adaptive recovery with exponential convergence rate 
recorded in Figure |H| An analogous filtering procedure can be carried out in the 
dual Fourier space, and as before, it hinges on a filter with an adaptive degree p 




Figure 5: Adaptive reconstruction of 
the piecewise smooth data from its us- 
ing N = 128-modes using (|4.1|l . 



Figure 6: Log error for an adap- 
tive spectral mollifier based on N ~ 
32, 64, and 128 modes. 
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Equipped with this toolkit to process spectral projections of piecewise smooth 
data, we turn to consider their time evolution. Godunov-type methods are based 
evolution of cell averages. Cell averaging is dissipative, but projections onto higher- 
order moments are not. The following example, taken from ^B] shows what go 
wrong with global projections. We consider the Fourier method for 27r-periodic 
inviscid Burgers' equation, dtV^(x, t)+d x SN{v%(-, t))/2 = 0. Orthogonality implies 
that ||«jv(-j is conserved in time, and in particular, that the Fourier method 

admits a weak limit, vpf(-,t) — 1 v(t). At the same time, v(t) is not a strong limit, 
for otherwise it will contradict the strict entropy dissipation associated with shock 
discontinuities. Lack of strong convergence indicates the persistence of spurious 
dispersive oscillations, which is due to lack of entropy dissipation. With this in mind, 
we turn to discuss the Spectral Viscosity (SV) method, as a framework to stabilize 
the evolution of global projections without sacrificing their spectral accuracy. To 
this end one augments the usual Galerkin procedure with high frequency viscosity 
regularization 

^-v N (x,t)+V N V x ■ f(v N (-,t)) =-N £ <r(^)«*(t)0*(aO, (4.2) 

|fe|<iV 

where cr(-) is a low-pass filter satisfying <j(£) > (|£| 2s — jf\ . Observe that the 
SV is only activated on the highest portion of the spectrum, with wavenumbers 
|fc| > 7 iV (2s - 1 )/ 2s . Thus, the SV method can be viewed as a compromise between 
the stable viscosity approximations corresponding to s — but restricted to first 
order, and the spectrally accurate yet unstable spectral method corresponding to 
s = 00. 



m(x> 




-1 -0.5 O 0.5 1 

Figure 7: Enhanced detection of edges 
with v(x) given by v(x) — —sgnx ■ 
cos(x + x ■ sgnx /2)lr_ 7] . ;7r i (x). 




-1 -0.5 0.5 1 



Figure 8: Legendre SV solution of in- 
viscid Burgers' equation. Reconstruc- 
tion in the direction smoothness. 



The additional SV on the right of (|4.2|) is small enough to retain the formal spectral 
accuracy of the underlying spectral approximation. At the same time, SV is large 
enough to enforce a sufficient amount of entropy dissipation and hence, to prevent 
the unstable spurious Gibbs' oscillations. The original approach was introduced 
in |18[ in conjunction with second-degree dissipation, s = 1. Extensions to several 
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space variables, non-periodic expansions, further developments of hyper SV methods 
with < s < oo, and various applications can be found in ^7]. We conclude with 
an implementation of adaptive SV method for simple inviscid Burgers' in Figure |H1 
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